library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.5.1     ✔ tibble    3.2.1
## ✔ lubridate 1.9.3     ✔ tidyr     1.3.1
## ✔ purrr     1.0.2     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
dados <- read_rds("../data/geomorfologia.rds")
glimpse(dados)
## Rows: 106
## Columns: 22
## $ sup     <chr> "I", "I", "I", "I", "I", "I", "I", "I", "I", "I", "I", "I", "I…
## $ solo    <chr> "LV", "LV", "LV", "LV", "LV", "LV", "LV", "LV", "LV", "LV", "L…
## $ amostra <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18,…
## $ x       <dbl> 0, 25, 50, 75, 100, 125, 150, 175, 200, 225, 250, 275, 300, 32…
## $ amg     <dbl> 0.2, 0.1, 0.7, 0.4, 0.4, 0.4, 1.2, 0.8, 1.1, 1.2, 0.1, 0.2, 0.…
## $ ag      <dbl> 3.72, 4.27, 5.00, 3.80, 3.10, 3.80, 3.60, 4.70, 4.50, 5.90, 5.…
## $ am      <dbl> 20.4, 22.6, 22.7, 23.7, 22.3, 23.8, 23.1, 25.8, 25.5, 32.8, 33…
## $ af      <dbl> 22.9, 23.6, 22.2, 24.4, 24.6, 19.1, 21.7, 21.1, 18.9, 19.8, 19…
## $ amf     <dbl> 30.0, 28.4, 26.9, 26.7, 26.9, 27.1, 26.5, 24.7, 25.4, 21.7, 19…
## $ silte   <dbl> 1.2, 1.2, 1.2, 0.6, 2.1, 2.2, 0.7, 0.2, 2.5, 0.2, 2.5, 2.6, 0.…
## $ argila  <dbl> 21.5, 20.4, 21.4, 20.5, 20.7, 23.5, 23.1, 22.7, 22.0, 18.5, 20…
## $ s_a     <dbl> 0.05, 0.05, 0.05, 0.02, 0.10, 0.09, 0.03, 0.01, 0.11, 0.01, 0.…
## $ af_ag   <dbl> 6.16, 5.53, 4.44, 6.42, 7.94, 5.03, 6.03, 4.49, 4.20, 3.36, 3.…
## $ p       <dbl> 42, 22, 41, 27, 11, 12, 11, 16, 38, 25, 25, 6, 6, 7, 5, 4, 4, …
## $ p_h     <dbl> 4.2, 3.8, 4.8, 4.0, 4.4, 4.0, 4.8, 5.4, 4.4, 5.2, 4.5, 5.1, 4.…
## $ k       <dbl> 0.27, 0.11, 0.34, 0.13, 0.11, 0.14, 0.23, 0.28, 0.19, 0.14, 0.…
## $ ca      <dbl> 1.4, 0.4, 2.4, 0.7, 1.4, 0.6, 1.6, 3.3, 1.6, 2.9, 1.3, 1.6, 0.…
## $ mg      <dbl> 0.3, 0.1, 0.4, 0.1, 0.3, 0.1, 0.7, 1.3, 0.5, 1.7, 0.6, 0.8, 0.…
## $ h_al    <dbl> 5.2, 5.8, 4.2, 5.2, 4.2, 5.2, 3.4, 2.5, 5.2, 3.1, 4.2, 2.5, 4.…
## $ sb      <dbl> 1.97, 0.61, 3.14, 0.93, 1.81, 0.84, 2.53, 4.88, 2.29, 4.74, 1.…
## $ t       <dbl> 7.17, 6.41, 7.34, 6.13, 6.01, 6.04, 5.93, 7.38, 7.49, 7.84, 6.…
## $ v       <dbl> 27, 10, 43, 15, 30, 14, 43, 66, 31, 60, 32, 50, 22, 35, 36, 34…
dados |>
  ggplot(
    aes(
      x=x,
      y=argila,
      color=sup
    ))+
      geom_point()

#Aplicando um moelo de delineamento inteiramente casualizado

y<-dados |> pull(argila)
trat<-dados |> pull(sup) |>as_factor()
trat
##   [1] I   I   I   I   I   I   I   I   I   I   I   I   I   I   I   I   I   II 
##  [19] II  II  II  II  II  II  II  II  II  II  II  II  II  II  II  II  II  II 
##  [37] II  II  II  II  II  II  II  II  II  II  II  II  II  II  II  II  II  II 
##  [55] II  II  II  II  II  II  II  II  II  II  II  II  II  II  II  II  II  II 
##  [73] II  II  II  II  II  II  II  III III III III III III III III III III III
##  [91] III III III III III III III III III III III III III III III III
## Levels: I II III

##Criar modelo para ánalise

mod<-aov(y~trat)
mod
## Call:
##    aov(formula = y ~ trat)
## 
## Terms:
##                     trat Residuals
## Sum of Squares  1255.866  2727.294
## Deg. of Freedom        2       103
## 
## Residual standard error: 5.145734
## Estimated effects may be unbalanced

##Estrutura do objeto mod

str(mod)
## List of 13
##  $ coefficients : Named num [1:3] 21.1 -6.34 -10.96
##   ..- attr(*, "names")= chr [1:3] "(Intercept)" "tratII" "tratIII"
##  $ residuals    : Named num [1:106] 0.4 -0.7 0.3 -0.6 -0.4 ...
##   ..- attr(*, "names")= chr [1:106] "1" "2" "3" "4" ...
##  $ effects      : Named num [1:106] -150.336 1.948 -35.385 -0.611 -0.411 ...
##   ..- attr(*, "names")= chr [1:106] "(Intercept)" "tratII" "tratIII" "" ...
##  $ rank         : int 3
##  $ fitted.values: Named num [1:106] 21.1 21.1 21.1 21.1 21.1 ...
##   ..- attr(*, "names")= chr [1:106] "1" "2" "3" "4" ...
##  $ assign       : int [1:3] 0 1 1
##  $ qr           :List of 5
##   ..$ qr   : num [1:106, 1:3] -10.2956 0.0971 0.0971 0.0971 0.0971 ...
##   .. ..- attr(*, "dimnames")=List of 2
##   .. .. ..$ : chr [1:106] "1" "2" "3" "4" ...
##   .. .. ..$ : chr [1:3] "(Intercept)" "tratII" "tratIII"
##   .. ..- attr(*, "assign")= int [1:3] 0 1 1
##   .. ..- attr(*, "contrasts")=List of 1
##   .. .. ..$ trat: chr "contr.treatment"
##   ..$ qraux: num [1:3] 1.1 1.11 1.16
##   ..$ pivot: int [1:3] 1 2 3
##   ..$ tol  : num 1e-07
##   ..$ rank : int 3
##   ..- attr(*, "class")= chr "qr"
##  $ df.residual  : int 103
##  $ contrasts    :List of 1
##   ..$ trat: chr "contr.treatment"
##  $ xlevels      :List of 1
##   ..$ trat: chr [1:3] "I" "II" "III"
##  $ call         : language aov(formula = y ~ trat)
##  $ terms        :Classes 'terms', 'formula'  language y ~ trat
##   .. ..- attr(*, "variables")= language list(y, trat)
##   .. ..- attr(*, "factors")= int [1:2, 1] 0 1
##   .. .. ..- attr(*, "dimnames")=List of 2
##   .. .. .. ..$ : chr [1:2] "y" "trat"
##   .. .. .. ..$ : chr "trat"
##   .. ..- attr(*, "term.labels")= chr "trat"
##   .. ..- attr(*, "order")= int 1
##   .. ..- attr(*, "intercept")= int 1
##   .. ..- attr(*, "response")= int 1
##   .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv> 
##   .. ..- attr(*, "predvars")= language list(y, trat)
##   .. ..- attr(*, "dataClasses")= Named chr [1:2] "numeric" "factor"
##   .. .. ..- attr(*, "names")= chr [1:2] "y" "trat"
##  $ model        :'data.frame':   106 obs. of  2 variables:
##   ..$ y   : num [1:106] 21.5 20.4 21.4 20.5 20.7 23.5 23.1 22.7 22 18.5 ...
##   ..$ trat: Factor w/ 3 levels "I","II","III": 1 1 1 1 1 1 1 1 1 1 ...
##   ..- attr(*, "terms")=Classes 'terms', 'formula'  language y ~ trat
##   .. .. ..- attr(*, "variables")= language list(y, trat)
##   .. .. ..- attr(*, "factors")= int [1:2, 1] 0 1
##   .. .. .. ..- attr(*, "dimnames")=List of 2
##   .. .. .. .. ..$ : chr [1:2] "y" "trat"
##   .. .. .. .. ..$ : chr "trat"
##   .. .. ..- attr(*, "term.labels")= chr "trat"
##   .. .. ..- attr(*, "order")= int 1
##   .. .. ..- attr(*, "intercept")= int 1
##   .. .. ..- attr(*, "response")= int 1
##   .. .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv> 
##   .. .. ..- attr(*, "predvars")= language list(y, trat)
##   .. .. ..- attr(*, "dataClasses")= Named chr [1:2] "numeric" "factor"
##   .. .. .. ..- attr(*, "names")= chr [1:2] "y" "trat"
##  - attr(*, "class")= chr [1:2] "aov" "lm"

##Extraindo os residuos do modelo

rs<-rstudent(mod)
rs
##            1            2            3            4            5            6 
##  0.079739246 -0.139552651  0.059803619 -0.119613529 -0.079739246  0.478958251 
##            7            8            9           10           11           12 
##  0.398994804  0.319106208  0.179436025 -0.518972767 -0.119613529 -0.379016227 
##           13           14           15           16           17           18 
## -0.199380742  0.458959546 -0.019934229 -0.099675805 -0.299143988  1.425307464 
##           19           20           21           22           23           24 
##  1.225272340  1.125799089  0.730837666  2.140207705  2.370928509  2.370928509 
##           25           26           27           28           29           30 
##  1.305103402  1.912650126  1.165548006  1.305103402  1.425307464  0.750492052 
##           31           32           33           34           35           36 
##  1.205350446  1.425307464  1.125799089  0.066039221  0.848895796  0.534710709 
##           37           38           39           40           41           42 
##  1.185442416  0.554292664  1.205350446  0.691553436  0.730837666  0.750492052 
##           43           44           45           46           47           48 
##  1.974419977  0.789826504  0.202556676 -0.148443870 -0.813952274 -0.499983185 
##           49           50           51           52           53           54 
##  0.007547180 -0.031446728  0.495564613 -0.519554627 -0.715627048 -1.110426371 
##           55           56           57           58           59           60 
## -1.389703543 -1.249711993 -1.429842586 -1.269665614 -0.774596123 -1.289633882 
##           61           62           63           64           65           66 
## -1.651851288 -0.774596123 -1.530485712 -1.570866517 -1.733147594 -0.912507996 
##           67           68           69           70           71           72 
## -1.189936578 -1.031119648 -0.187453873 -1.170039114 -0.813952274 -0.480417361 
##           73           74           75           76           77           78 
## -0.480417361 -1.409764874 -1.530485712 -1.209847718 -1.490176533 -1.189936578 
##           79           80           81           82           83           84 
## -1.329615342 -0.403227845 -0.265046270 -0.978988707 -0.324244221  0.089781704 
##           85           86           87           88           89           90 
## -0.919081451 -0.978988707 -1.159352014 -0.978988707 -0.939039543 -0.879197047 
##           91           92           93           94           95           96 
## -0.008758848 -0.186157061 -0.205875607  1.895996678  0.464709923  0.267238263 
##           97           98           99          100          101          102 
##  0.267238263  0.981209424  1.444431251  0.188347901  0.444942530  0.781869685 
##          103          104          105          106 
## -0.700194582  0.484482647  0.801759078  0.841565924

##Extraindo os preditos pelo modelo

yp<-predict(mod)
yp
##        1        2        3        4        5        6        7        8 
## 21.10000 21.10000 21.10000 21.10000 21.10000 21.10000 21.10000 21.10000 
##        9       10       11       12       13       14       15       16 
## 21.10000 21.10000 21.10000 21.10000 21.10000 21.10000 21.10000 21.10000 
##       17       18       19       20       21       22       23       24 
## 21.10000 14.76129 14.76129 14.76129 14.76129 14.76129 14.76129 14.76129 
##       25       26       27       28       29       30       31       32 
## 14.76129 14.76129 14.76129 14.76129 14.76129 14.76129 14.76129 14.76129 
##       33       34       35       36       37       38       39       40 
## 14.76129 14.76129 14.76129 14.76129 14.76129 14.76129 14.76129 14.76129 
##       41       42       43       44       45       46       47       48 
## 14.76129 14.76129 14.76129 14.76129 14.76129 14.76129 14.76129 14.76129 
##       49       50       51       52       53       54       55       56 
## 14.76129 14.76129 14.76129 14.76129 14.76129 14.76129 14.76129 14.76129 
##       57       58       59       60       61       62       63       64 
## 14.76129 14.76129 14.76129 14.76129 14.76129 14.76129 14.76129 14.76129 
##       65       66       67       68       69       70       71       72 
## 14.76129 14.76129 14.76129 14.76129 14.76129 14.76129 14.76129 14.76129 
##       73       74       75       76       77       78       79       80 
## 14.76129 14.76129 14.76129 14.76129 14.76129 14.76129 14.76129 10.14444 
##       81       82       83       84       85       86       87       88 
## 10.14444 10.14444 10.14444 10.14444 10.14444 10.14444 10.14444 10.14444 
##       89       90       91       92       93       94       95       96 
## 10.14444 10.14444 10.14444 10.14444 10.14444 10.14444 10.14444 10.14444 
##       97       98       99      100      101      102      103      104 
## 10.14444 10.14444 10.14444 10.14444 10.14444 10.14444 10.14444 10.14444 
##      105      106 
## 10.14444 10.14444

Começar o diagnostico construindo um arquivo cm sup,y,yp,rs

diagnostico<-tibble(
  trat,
  y,
  yp,
  rs
)
diagnostico
## # A tibble: 106 × 4
##    trat      y    yp      rs
##    <fct> <dbl> <dbl>   <dbl>
##  1 I      21.5  21.1  0.0797
##  2 I      20.4  21.1 -0.140 
##  3 I      21.4  21.1  0.0598
##  4 I      20.5  21.1 -0.120 
##  5 I      20.7  21.1 -0.0797
##  6 I      23.5  21.1  0.479 
##  7 I      23.1  21.1  0.399 
##  8 I      22.7  21.1  0.319 
##  9 I      22    21.1  0.179 
## 10 I      18.5  21.1 -0.519 
## # ℹ 96 more rows
#Normalidade dos residuos

diagnostico |> 
  ggplot(
    aes(
      x=rs
    )
  ) +
  geom_histogram(
    bins=10,
    color="black",
    fill="gray30"
    )+
      theme_bw()

#QQplot
diagnostico |> 
  ggplot(
    aes(
      sample=rs
    )
  )+
  stat_qq()+
  stat_qq_line(
    color="blue"
  )

diagnostico |>
  ggplot(aes(sample = rs)) +
  stat_qq() +
  stat_qq_line(color = "blue") +
  facet_wrap(~ trat) +
  theme_bw()

library(plotly)
## 
## Anexando pacote: 'plotly'
## O seguinte objeto é mascarado por 'package:ggplot2':
## 
##     last_plot
## O seguinte objeto é mascarado por 'package:stats':
## 
##     filter
## O seguinte objeto é mascarado por 'package:graphics':
## 
##     layout
# Calcular quantis teóricos e empíricos
qq_data <- qqnorm(diagnostico$rs, plot.it = FALSE)
# Suponha que 'trat' seja um fator indicando grupos
# Converta em numérico apenas para ilustrar no 3D
group_numeric <- as.numeric(diagnostico$trat)

plot_ly(
  x = ~qq_data$x,
  y = ~qq_data$y,
  z = ~group_numeric,
  type = 'scatter3d',
  mode = 'markers'
) %>%
  layout(
    scene = list(
      xaxis = list(title = "Quantis Teóricos"),
      yaxis = list(title = "Quantis Empíricos"),
      zaxis = list(title = "Grupo")
    )
  )
# install.packages("nortest")
#Aplicando os testes de normalidade
shapiro.test(rs)#nao rejeitamos H0
## 
##  Shapiro-Wilk normality test
## 
## data:  rs
## W = 0.97323, p-value = 0.03052
nortest::lillie.test(rs)#nao rejeitamos H0
## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  rs
## D = 0.072002, p-value = 0.1941
nortest::cvm.test(rs)#nao rejeitamos H0
## 
##  Cramer-von Mises normality test
## 
## data:  rs
## W = 0.077302, p-value = 0.2222
# Conclusão
## Os resíduos seguem uma distribuição NORMAL
##Gráfico dos 5 números ou Boxplot

diagnostico |> 
  ggplot(aes(
    x=trat,
    y=y,
    fill=trat,
  )
  )+
    geom_boxplot()

diagnostico |>
  ggplot(aes(x = trat, y = y, fill = trat)) +
  # Gráfico de violino com ajuste na largura
  geom_violin(trim = FALSE, alpha = 0.7, scale = "width") +
  
  # Adicionar pontos individuais com jitter para visualizar a distribuição
  geom_jitter(width = 0.1, alpha = 0.5, size = 1.5, color = "black") +
  
  # Adicionar a mediana de forma destacada
  stat_summary(fun = median, geom = "crossbar", 
               width = 0.3, color = "white", 
               fatten = 2, linetype = "solid") +
  
  # Ajuste de escala de cor (ex: viridis - requer library(viridis) ou scale_fill_viridis_c)
  scale_fill_brewer(palette = "Pastel1") +
  
  # Título e rótulos
  labs(title = "Distribuição por Tratamento", 
       x = "Tratamento", 
       y = "Resposta",
       fill = "Tratamento") +
  
  # Tema minimalista
  theme_minimal(base_size = 14) +
  theme(
    legend.position = "top",
    plot.title = element_text(hjust = 0.5, face = "bold")
  )

#Teste de homocedasticidade
#Teste de Levene
lawstat::levene.test(y,trat)
## 
##  Modified robust Brown-Forsythe Levene-type test based on the absolute
##  deviations from the median
## 
## data:  y
## Test Statistic = 22.666, p-value = 6.96e-09
#Bartlett
bartlett.test(y,trat)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  y and trat
## Bartlett's K-squared = 31.396, df = 2, p-value = 1.522e-07
##estudar a regularidade
##calcular media e variancia dos valores de y em funçao dos trat
diagnostico |> 
  group_by(trat) |> 
  summarise(
    log_media=log(mean(y)),
    log_variancia=log(var(y))
  ) |> 
  ggplot(
    aes(
      x=log_media,
      y=log_variancia
      )
    )+
  geom_point()+
  geom_smooth(
    method = "lm",
              se=FALSE
    )
## `geom_smooth()` using formula = 'y ~ x'

#Análise de regressão
modelo_linear<-lm(log_variancia~log_media,
                  data=diagnostico |> 
                    group_by(trat) |> 
                    summarise(
                      log_media=log(mean(y)),
                      log_variancia=log(var(y))
                      )
)
#transformação Box and Cox